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ABSTRACT 


N A Computer Model of Solar Panel -Plasma Interactions" 


High-power solar arrays for satellite power systems are presently being 
planned with dimensions of kilometers, and with tens of kilovolts distributed 
over their surface. Such systems will face many plasma Interaction problems, 
such as power leakage to the plasma, particle focusing, and anomalous arcing 
to name a few. In most cases, these effects cannot be adequately modeled 
without detailed knowledge of the plasma-sheath structure and space charge 
effects. This report details the work performed under contract NAS 9-15796 to 
adapt the computer program PANEL to augment the laboratory studies of a 1 x 10 
meter solar array In a simulated low Earth orbit plasma being conducted In the 
chamber A facilities at NASA/Johnson Space Center. The plasma screening 
process is discussed, program theory is outlined, and a series of calibration 
models Is presented. These models are designed to demonstrate that PANEL Is 
capable of accurate self-consistent space-charge calculations. Such models 
Include PANEL predictions for the Child-Langmulr diode problem. Also Included 
are two models of the Interactions of an Infinitely long one meter wide solar 
array In a dense, 10 eV plasma. 
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1: INTRODUCTION 


The Interaction of a large high voltage solar array with a space or 
laboratory plasma cannot. In general, be modeled analytically. For this 
reason, a computer program, PANEL, has been developed to calculate potentials, 
densities, and currents on a three-dimensional grid of points. A major 
feature of PANEL Is that It Includes space-charge effects In a self-consistent 
manner. The method used In PANEL Is an extension to three dimensions of the 
Inside-out method developed by Parker (1964), and used by Parker and Whipple 
(1970) to model two-electrode probes on a satellite. More recently Parker 
(1976, 1977a) has used the method to calculate sheath and wake structures 
about disk and pill-box shaped objects In flowing plasmas. An early version 
of PANEL, written by Parker (1977) was used by Relff, Freeman, and Cooke 
(1980) to model the Interaction of a geosynchronous substomt plasma with the 
NASA/Marshall Space Flight Center baseline design for the solar power satel- 
lite (Hanley, 1978). The purpose for the further development of PANEL, 
reported here, has been to produce a code capable of augmenting the laboratory 
studies of a 10 meter solar array In a simulated low Earth orbit plasma being 
conducted at NASA/Johnson Space Center (McCoy and Konradl, 1978). 

Section 2 of this report contains a general discussion of the plasma 
screening process, applications of the Chlld-Langmulr diode law, and a 
development of an analytic model for comparison to PANEL results. Section 3 
describes the Inside-out method as It Is used In PANEL to solve the coupled 
Vlasov-Polsson equations In two- and three-dimensions. And In section 4, 
current PANEL results are presented. Appendix A contains a review of the 
symbols and units used In this report. Appendix 6 gives a subroutine linkage 
chart, and Appendix C outlines the modifications to PANEL for the two- 
dimensional option. The emphasis so far, has been on testing PANEL against the 
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analytic models discussed In section 2. Also Included are two two-dimensional 
models of an "Infinitely long" charged panel. 


2: THE PLASMA SHEATH 

Perhaps the best known example of plasma screening Is the Debye treatment 
of the plasma screening of an Isolated test charge. A poslti/e test 
charge, 6Q, placed In a plasma of temperature T, will attract electrons and 
repel Ions so as to develop a surrounding sheath with a potential distribution 
given by. 


»('•) * *FTr exp <- r/ V 

0 

where Xq * (e 0 kT/N 0 e 2 ) 1 / 2 Is the Debye length. Implicit in the derivation of 
this equation (Jackson, 1962) are the assumptions that the charge has 
negligible cross-section, and that V(r) « kT/e for r > Xg. For a microscopic 
body of radius R, satisfying these assumptions, we can write 

V R 

V(r) - exp ( *r/x D ) (2-1) 

where Vj, Is the surface potential of the body. 

For macroscopic bodies with a high degree of symmetry, sheath structures 
may be calculated by capitalizing on the constants of the motion allowed by 
the symmetry, e.g., angular momentum (Whipple, 1977). In general though, 
self-consistent treatment of a macroscopic body requires computer modeling. 
In spite of this difficulty, a better understanding of the shielding process 
can be gained by studying current limiting by space charge In the 1-D planar 
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electron diode. 

The first theoretical treatments of the electron diode were published 
Independently by Child (1911) and Langmuir (1913). Variations on this problem 
have been studied by Fay et al. (1938), and the general topic of space charge 
effects In vacuum tubes Is treated In the book by Blrdsall and Bridges (1966). 

Consider the three electrode system shown In Figure 2/1. At x ■ -d, we 

have a cathode, with zero potential capable of emitting unlimited quantities 
of electrons all with zero velocity. At x ■ 0, we have a transparent screen 
at potential V Q , and at x ■ Xj, we have a non-emitting anode at potential V . 
The kinetic energy of an electron at x Is 

■| mgV 2 - eV. (2-2) 

The current density for electrons at x Is 

J * pv (2-3) 

where p Is the charge density. Poisson's equation In one dimension Is 


d 2 V 
dx 2 " 


p/e 0 ‘ 


(2-4) 


Substituting for p from equation (2-3) and then for v from 
have 

d 2 V f m e)i/2 

dx 2 ‘ e o • 


equation (2-1), we 
(2-4) 


Multiply this by 2^- and Integrate from (0, V Q ) to (x, V x ); 

(f i ) 2 i; - - jj • (Vx‘ /J - V/ 2 ) (2-5) 

where b 2 ■ (2e/m e ) l/r2 - 2.336 k 10- 6 (amps/volts 3 / 2 ). The boundary condl- 

dV 

tlon that we desire at x ■ 0, Is E x ■ - * 0, a common definition of the 
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sheath edge In the spacecraft charging problem. Using this condition, and 
taking the positive square root of equation (2-5), 


£ ( 2 - 6 ) 
where J ■ -J. Solving *or dx we Integrate again from (0, V Q ) to (x , V ); 

e u 1 1 

x, . -jj= V V- [1 ♦ 2 (^) l/2 l . [1 - $*) l/2 ] 1/2 . (2-7) 

This equation can be applied to region I, -d < x < 0, where V Q - 0 (zero 
Initial velocity) to recover the Chlld-Langmulr (C-L) result 

d 2 - 2.336 x 10- 6 V 3/2 /J e , (2-8) 


where d 2 and J e must have the same unit of area. If d and V are fixed, 
equation (2-8) gives the maximum conducted current despite an unlimited supply 
of electrons. If d, V, and J e are all considered Independent, the sheath edge 
electric field that was set to zero will become the dependent variable. 

To apply equation (2-7) to a planar spacecraft surface at potential V, we 
identify Xj with the spacecraft surface, and x ■ 0 with the sheath edge where 
E * 0. Region I Is now Identified with the undisturbed plasma where V 0 repre- 
sents the average thermal energy of the electrons. Therefore, a transfor- 
mation of voltage Is needed because we want the sheath edge to be at zero 
potential; thus, * V ♦ V 0 


* 


V 

r 0 




define 


and 


C+ + D" 1 - 


(2-9) 


Substituting this last result Into equation (2-7) we get for the sheath 
thickness a, 

x ■ a ■ b (V 3/l 7/d£) • S(a) (2-10) 

s(*> ■ id ♦ • id ♦ 2 a + *r 1/2 ) • a - a + *r l/ 2 ) 1 / 2 ]i. (2-id 

Equation (2-10) Is written such that S(*) Is the correction to the usual 
Chlld-Langmulr result due to non-zero Initial velocities. Equation (2-11) Is 
plotted In Figure 2/2 as a function of log ( ♦" 1 ) • When applied to a 

Maxwellian plasma, S(*) will be only qualitatively correct since there will be 

a distribution of Initial velocities, but It should be reasonably accurate for 
larger values of ♦. 

Another variation of this problem Is given by Blrdsall (1966). The 
conditions are Illustrated In the lower portion of Figure 2/1, with the grids 
at x ■ 0 and X! both at the same positive potential V l# and the separation 
distance xi considered fixed. The negative space charge of the electrons In 
the gap between zero and xi will depress the potential In the gap and give 
rise to current limitation If the potential drops to zero. This variant Is 

more suited for comparison to PANEL, since the geometry Is fixed and only 

voltages and charge densities vary. The potential distribution In the gap Is 
determined by subdividing region II Into regions A and B whose boundary at x m 
is the point of minimum potential where we have also the condition of zero 
electric field. The potential can be obtained separately In each region with 
exactly the same approach that led to equation (2-7) to give 


D[(2 - f) J]'* /2 • (V A l/2 ♦ 2V m ‘ /2 ] * ( V A 1/2 - V ffl I/2 ) l/2 ■ % - x A 

(2-12) 

b(fj)* 1/2 • (V B l/2 ♦ 2V m 1/2 ) • (V B 1/2 - V m 1/2 ) l/2 ■ x B - x,, 

The factor, f, Is the fraction of transmitted current; If V m > 0, f ■ 1. 
Equations (2-12) may be solved for by setting (V A , x A ) ■ (V;, 0) and 
( V B» X B^ " (V i . X;). For the case V m > 0 and f ■ 1, we find x m ■ *i/2, and 
equations (2-12) can be combined to give Blrdsall's equation: 

(♦»/* - V/2) (* 1/2 ♦ 2«» l/2 ) 2 * 6(t - 1/2) 2 (2-13) 

where ♦ ■ V/Vj, ■ Vjn/Vg and 5 * x/xj. The dimensionless current 8, Is 

defined by 8 ■ Je/b 2 V; 3 / 2 X;" 2 , where the normalizing current Is C-L current 
for a diode with separation x , and potential drop V t . The value of the mini- 
mum potential for a given current Is found by evaluating equation (2-13) at 
x ■ 0, giving 


4 *m /2 ' - 1) ■ 0. (2-l«) 

Figure 2/3a (from Blrdsall) gives as a function of Input current 8. The 
dotted portion of the curve, has been labeled the "C-overlap" by Fay et al. 
(1938) and has been shown to contain no stable solutions (Blrdsall, 1966). 
Figure 2/3b (also from Blrdsall) gives selected potential profiles from 
equation (2-13) for the range 0 < 8 < 8. 

For solutions with 8 > 8 we set ■ 0 In equation (2-12). Allowing now 
for f < 1 and a non-symmetrlc potential distribution, we can derive Blrdsall's 
dimensionless equations; 


C 


(2-15) 


. 3/2 


•a ■ (2 - f) • (c - c m ) 


*b /2 ■ ft (c - w 2 . 


. rf (2 - fn 1/2 - f 
** 2(1 - f) • 

. . 2 r Cf (2 • f)l| /2 ♦ l i 

0 2 L f (2 - f) J 


(2-16) 

(2-17) 


where € m • Xg/xi* Figure 2/3c (also borrowed from Blrdsall), shows a few 
selected potential profiles from equations (2-15). 

The multiplicity of solutions In the region 4 < 6 < 8 means that there 
should be hysteresis In the behavior of the classical diode model. This 
hysteresis was observed experimentally by Gill (1925). He also observed the 
predicted current limiting. It should be mentioned that although the condi- 
tion V m ■ 0 does lead to mathematical solutions for the range 4 < S < • these 
solutions are a result of the enforced time- Independence of the classical 
method and can be shown to be unstable when time dependence Is considered 
(Blrdsall, 1966). This lack of stability has been confirmed by the experi- 
ments of Salzberg and Haeff (1938). 

These solutions of what I call the gap problem have been very useful In 
the development of PANEL. PAPCL's predictions for the gap model will be 
presented In section 4. 

Returning now to the spacecraft plasma sheath problem. It Is common to 
use the C-L sheath thickness given In equation (2-8) as an estimate of 
satellite sheath thickness for appropriate conditions. These conditions are: 


7 


1. Satellite dimensions should be much larger than all estimates of 
sheath thickness, such that a planar approximation Is justified. 

2. The surface potential Is nuch greater than the plasma temperature, so 
the Initial velocities of particles entering the sheath can be neglected 
[or accurately accounted for by equation (2-10)]. Also repelled 
particles must not penetrate significantly Into the sheath since the C-L 
treatment considers only attracted particles. 

3. The current Is assumed to be the random thermal current 
(Jq ■ N 0 e /kT/Zwm) of attracted particles falling on the sheath edge. 

For diode geometries other than planar, Langmuir (1913) has shown that 
the space charge limited current will always be proportional to V 3 ^ 2 , however, 
the distribution of potential In space does depend on geometry. The problems 
of current flow between concentric spheres and cylinders has been addressed by 
Langmuir and Blodgett (1924). Their solutions take the form of equation (2-8) 
with d replaced by various series expansions In terms of the ratios of the 
electrodes' radii, with the results presented In tabular form. Parker (1979) 
has adapted these results to estimate sheath thickness for charged spherical 
satellites, and provides a convenient fit to those results. In the following 
equations a ■ sheath rdlus, r Q ■ body radius, and d Is the C-L screening 
distance given by equation (2-10) with ♦ ■ •, or by equation (2-20) below: 
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It Is Interesting to stop at this point and compare screening length 
estimates. On one hand, we have the Debye length for the microscopic case 
with low potentials 


x d ' [0] I/2 "» t * r ‘ »-») 
and, on the other hand, the planar C-L screening length 

dd V s/4 • 9.34 (T(ev) ^(cm- 3 ))‘ l/2 V 3/4 meters (2-20) 

J Ng*kT 

where 1 have substituted for the current the thermal current J 0 . Note the 
difference In the temperature and density dependencies; (T/N) ' as opposed to 
(N*Tr l/ \ A brief example will help demonstrate the Inappropriateness of 
applying the Debye model to large objects. Consider an object of radius 
r 0 • 10m, at a potential of 100V In a plasma with T # ■ 1 ev, and N ■ 100/cc. 
For these conditions, the C-L distance Is by (2-20), d - 29.5 m, and from 
(2-18) we have a * 2.4 r Q ■ 24m; whereas with a Debye length of 0.74 m, 
equation (2-1) predicts a ■ 13 m for a sheath edge potential of one volt. The 
Debye model predicts screening on the order of a few Debye lengths which 
significantly underestimates the sheath thickness. 

The planar, one-dimensional C-L screening length has Its own short- 
comings. In applying It to space conditions one assumes that the screening 
length will be small compared to surface dimensions so that the surface can be 
approximated as Infinite. But this Is Incompatible with the assumed boundary 
condition of an undisturbed plasma one screening length away from the surface. 
Since repelled particles will be reflected at the sheath edge, the repelled 
particle distribution will be nearly Isotropic. On the other hand, the 
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attracted parti cits coma to tha edge from one direction only, resulting In 
only a hemispherical distribution out to a point where the surface no longer 
appejrs Infinite. Such a plasma cannot even be quasi-neutral unless the 
potential at the "sheath edge" Is significantly »,:>n-zero, which would be 
contrary to the original conditions of the C-L model. This difficulty Is not 
as severe for the cylindrical and spherical extensions. This situation can be 
Improved somewhat If we relax the definition of the sheath edge to require 
only that the electric field be near zero there, and let the potential deviate 
from zero (closer to the surface potential) so as to draw In more attracted 
particles and -educe the repelled particle density. This In turr requires 
that wi Invoke a presheath region to match the sheath edge to regions where 
the attracted and repelled distributions are Identical, and the potentl'1 and 
electric fields vanish. 

This presheath problem was recognized and accounted for by Langmuir 
(1929) In his analysis of a plasma discharge between plane parallel 
electrodes, however the resolution of that problem Is too Involved for 
presentation here, and net directly applicable to the problem of a planar 
satellite In a colllslonless plasma. The presheath problem has also been 
studied by Parker (1980) for an extremely large spherical body (large compared 
to all estimates of sheath thickness) In a colllslonless plasma, using a 
technique similar to the Inside-out method but solving a condition of 
quasi neutrality Instead of Poisson's equation. The sheath was modeled as a 
potential discontinuity at the surface. Results Indicate that no matter how 
thin the sheath gets, the presheath thickness will always be comparable to 
body dimensions (as one would expect from geometrical shadowing conside- 
rations). In the presneath region the potential drop Is of order kT/e and the 


Ion and electron densities are essentially equal, but reduced from ambient 
values. One-dimensional theory predicts that such a body will collect a 
current density of attracted particles equal to the random thermal density 
current, but a conclusion of Parker's study was that the collected current 
density will be Increased in the ~*esheath region by a factor dependent upon 
the body potential but approaching a limiting value of 1.45 for Infinite body 
potential, Independent of body shape. 

Another complication Is presented by secondary and photoelectron 
emission. The extent to which these additional sources will modify a plasma 
sheath is of course dependent upon the emission flux. Guernsey and Fu (1970) 
and Fu (1971), have studied the case charactarized by (Ny/Ng) > (T e /Ty) > 1 
where the u subscript refers to the photo or secondary electrons, and the e 
for the plasma electrons. We would expect these conditions to lead to a 
positive surface potential of a few volts and a monotonic decrease to zero 
with increasing distance from the surface. Their study confirmed that 
solution, but also revealed the existence of another non-monotonlc type with a 
negative "overshoot" potential minimum located about one photoelectron Debye 
length from the surface. This overshoot was also accompanied by a lowering of 
the equilibrium surface potential by an amount roughly equal to the overshoot 
potential. Without a time dependent analysis one cannot decide which solution 
is the true steady-state, but energy considerations suggested that the non- 
monotonic solution is the true steady-state when (^y) * where H y is 
the peak in the photoelectron kinetic energy spectrum. 

If the preceeding picture of the plasma screening process seems incom- 
plete, it is for a few good reasons. One is that many aspects of the problem 
have not been fully investigated yet, but more importantly, there are just too 
many parameters to ever hope to develop a "one-size fits air theory of plasma 
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screening. In n\y opinion, the objectives of this line of research Is to 
develop the necessary methods and computational tools to solve specific 
problems as they arise. 




3: THE INSIDE-OUT METHOD 

The classical theory of electrodynamics states that the scalar 
electrostatic potential V(x) and the charge density p(V(x)) will satisfy 
Poisson's equation 

7 2 V(J) - -p(V(J))/c 0 . (3-1) 

In problems where the charge density does not depend upon the potential, 
equation (3-1) becomes an inhomogeneous linear elliptic partial differential 
equation. For such equations, the theory of partial differential equations 
(Jackson, 1962), will guarantee a unique solution Interior to a closed 
boundary S, on which is specified either (but not both) the potential V(x $ ) 

A 

(Dlrlchlet boundary conditions), or the normal derivative 3V(x $ )/3n s (Neuman 
boundary conditions). Unique solutions may also be obtained for problems with 
mixed boundary conditions with Dlrlchlet conditions on part of the boundary, 
and Neuman for the rest. For the general non-linear problem where the charge 
density depends on the distribution of potential, there are no uniqueness or 
existence guarantees for solutions to equation (3-1). Experience, however, 
leads us to believe that the physically real problems that we encounter In the 
study of plasma screening do have at least one self-consistent solution for V 
and p. It Is this experience that leads us to pursue solutions to such 
problems. 
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The Inside-out method adopts an Iterative approach to solving plasma 
sheath problems. The best estimate for p(x) Is used In equation (3-1) to 
obtain a new estimate for V(x). Next, new estimates for o(x) are obtained via 
the Vlasov equation and the latest values for V(x) and the process repeated. 
The calculation of p(x) has been labeled the "Vlasov problem" and the problem 
of finding solutions to (3-1) Is called the “Poisson problem". 

PANEL has the feature of being able to operate In both a two-dimensional 
and a three-dimensional mode. The three-dimensional version is presented 
first, and the conversions to two-dimensional operation are found In Appendix 
C. 


THE POISSON PROBLEM 

With p(x) temporarily considered known and Independent of V(x), equation 
(3-1) becomes linear, thus a well posed boundary value problem will have a 
unique solution. PANEL uses a standard finite-difference method to solve 
Poisson's equation (Collatz, 1960). The approach is to discretize the space 
to be modeled by constructing a three-dimensional grid of points j ic* An 
x-y plane at constant z is Illustrated In figure 3/1. The standard approach 
is to let the x, y, and z spacings all be a constant h so that there is a cube 
of volume h J associated with each Interior point. But, in modeling many 
objects it Is convenient to use variable spacing to achieve greater economy by 
allowing a higher density of points where a need is anticipated. This means 
that each Interval must be calculated, but the symbol h will still be used to 
represent a typical Interval. With variable spacing, the volume associated 
with a point P becomes a rectangular parallelepiped with faces located at the 
midpoints between P and its neighbors. The shaded area in figure 3/1 
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represents the x-y projection of this volume. Also Indicated in the figure Is 
the sense of the directions represented by the notation, N, S, E, W, U, D for 
north, south, east, west, up, and down, respectively. 

On this grid, we now develop a difference equation to aproxlmate (3-1). 
We start with the central difference operator, «, which Is defined as 

6 f (x) - f (x + h/2) - f (x - h/2). 

Applying this operator twice we get, 

« 2 f (x) » 6 [f (x + h/2) - f (x - h/2)] 

- [f (x ♦ h) - f (x) ] - [f(x) - f(x - h)] (3-2) 

* f(x + h) - 2f(x) + f(x - h). 

The connection between this last expression and the second derivative can 
be observed by first writing the Taylor series for f (x + h), 

f (x i h) - f(x) i » 37 (*) + ?r 0 (X) i yr h 2 0 (x) ♦ (3-3) 

substituting these series into (3-2), we find 


« 2 f(x) 


f(x h) -2f(x) + f(x - h) 

2 t h2 ir0 (x) + ir 1,4 0 (x) l + °< h6 >> 


(3-4) 


d 2 f 

and solving for — - 
dx z 
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(3-5) 


dx 2 h 2 « dx** 


The simplest approximation that we can make for the second partial derivative 
Is thus. 


a 2 f(x.y,z) a f(x » h, y,z) - 2f(x,y,z) » f(x - h,y,z) 
ax 2 h 2 


(3-6) 


with an error of order h 2 . Higher order approximations can be obtained by a 
similar Taylor series analysis (Collatz, 1960), but the resulting formulas 
sufficiently complicate the consideration of boundary conditions enough to 
discourage their use in favor of just reducing h as much as possible. This of 
course Increases the number of points required to model a given object; so If 
a machine size limit is reached and greater accuracy Is still required, the 
use of higher order approximations would be an option. 

To investigate the effect of using variable spacing we can let h ♦ h + , h. 
where both are positive numbers. Starting again at (3-3) with this change, we 
see that the odd order contributions to (3-4) no longer cancel, thus (3-5) 
becomes 


d 2 f ZlZf z ( h + - h -) df . , 

dx 2 ' (h+ 2 ♦ h. 2 ) * (h+ 2 ♦ h. 2 ) <1* ' ’ 


(3-7) 


so, if h + and h_ are not nearly equal, accuracy will be reduced. 

We could now formally construct a differenced form of Poisson's equation 
from (3-6), however It is more honest to present PANEL'S Poisson algorithm as 
originally developed by Parker (19/7b). We first throw Poisson's equation 
into partially dimensionless form by dividing by kT/e, so with • ■ Ve/kT and 
x 0 2 * e 0 kT/N 0 e 2 we get 
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»**(x) ■ X D **(n, • ih) ■ R, 


(3-8) 


where n e and are the electron and Ion densities In units of the ambient 
density Nq. Integrate now (3-8) over the cell volume associated with point P, 
and apply the divergence theorem to the left hand side; 

J/Jv 2 e d 3 x - f s ds - /// R d 3 x - Q. (3-9) 

where ^ Is the outward normal derivative at the surface of the cell. Q can 
be Identified as the net charge within the cell, however, this Identification 
Is not Implicit In the formal development. We next approximate the surface 
Integral In (3-9) by the sum: 


I Ap(-^)p s Q, (3-10) 

where F ■ N, S, E, W, U, D, and Ap Is the area on each of these faces. These 
areas are given by, 


a n * A s -T< x in - - Z k-J 

a e (WW 

*U ,A D 'T tx 1*l - X 1-l> ^>1 -yj-l)- 
The partlals (-j^J are approximated by the difference quotients: 


1*1 - *" ~ * fiii - *« ' 

3 " j n y J+l - w s yj-yj. 


(3-11) 


(3-12) 


and similarly for the E, W, U, and 0 directions, where ♦ Is the potential at 
the point P, and * N> * s , etc. are the neighboring potentials. Thus substi- 
tuting equations (3-11) and (3-12) Into (3-10) we obtain the algebraic 
expression, 

C N*N + C S # S + ^E*E + + Sj*U * S # D “ " Q» (3-13) 

where c . ^ x 1*l ' x 1-l^ z k*i ' Vi 1 

N 4(y 1+1 - yi) 


and likewise for tLrough Cq; C » l Cp. 

Equation (3-13) can be applied to each Interior point In the model, but 
exterior or boundary points require a modified treatment so as to Include the 
required boundary conditions. The types of boundary conditions (B.C.) used In 
PANEL are: 

1. Floating, where the outward normal derivative on the cell and model 
boundary is linearly related to the potential on the boundary*. 


*In the theory of boundary value problems. Independent specification of 
the normal derlvat* e and potential Is an over specification of the boundary 
conditions and there will be no solution unless the solution was already known 
and used to specify the B.C. Here we are specifying only a relation between 
the two conditions, but even this implies a knowledge of the Green function 
for the problem. For the case where the boundary Is far enough away from the 
"object" for the object to look like a point charge or at least a uniformly 
charged sphere, we can assume a Green function of 1/r, so we have the 
relations: 


* n • 
dn 


-n 


♦ . 


For a closer boundary, the possibility exists for finding the appropriate 
Green function, but this has not been pursued far enough to produce a useful 
algorithm. 
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2. Neuman, where the Inward normal component of the electric field Is 
specified. 

3. Olrlchlet, where the boundary potential Is specified. 

4. Extended Olrlchlet, where a boundary potential of zero Is assumed to 
exist one Interval beyond the usual model boundary. 

5. Reflection, where like condition 4, an extended boundary Is assumed, 
but with a potential equal to the nearest Interior neighbor. 

When a boundary Is assumed to represent "Infinity", l.e. a source of 
undisturbed plasma at zero potential, the boundary should be far enough away 
from the "object" that all boundary conditions give the same results. It Is 
frequently not possible to make grids that large so It becomes necessary to 
choose the B.C. which best approximates "Infinity" on a limited grid. Parker 
and Sullivan (1969) has addressed this problem, and concluded that for the 
spherical diode problem, the floating B.C. (1) produced the best "Infinity" 
approximation with the least computing time. The zero gradient B.C. (2) 
produced an effective infinity at a distance comparable to the floating B.C., 
but required about twice the computing time as B.C. (1). The zero potential 
B.C. (3) required a more distant boundary to produce similar results, and the 
required computing time was between that required for B.C. (1) and B.C. (2). 

All of these boundary conditions are effected by treating a boundary 
point as an Interior point, and by adding the appropriate "off-grid" poten- 
tial. Equation (3-12) and (3-13) also require modification at exterior points 
since If the maximum (minimum) value of x In the model Is Xjj ( Xj ) the point 
X II+i ^ X o ) does not The required modifications for A N , Aj, (3*/3n)^ 

and (aa/3n] w are respectively: 


A N * 4 ( x 1 + i * x 1-l^ z k+l ' Z k-J * 2 ( X II ’ x II-i^ z k+i ' z k-i^ for 1 


II; 
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A S ■ 7(*1 ■ x l^ I k+l ' z k- for 1 ■ l! 
X 1+ l - x< -(Xj Nl - x IX ) • 

(il) - V* ,*W, 

'in'n *1 - *1-1 1*2 - *iJ 


where the asterisk Indicates that •* Is chosen subject to the boundary 
condition. 

With the appropriate consideration of exterior points, we can now apply 
equation (3-13) to all grid points giving a system of linear equations that Is 
solved by the method of over-relaxation (O.R.) (Stlefel, 1963). Faster and 
more sophisticated methods are discussed by Hockney (1965), but O.R. has been 
chosen for Its programming simplicity and versltlllty. To derive the O.R. 
formula used In Panel's relaxation algorithm, we first cast the system of 
equations produced by equation (3-13) Into the form 

M 

I C ♦ - Q « 0, p « 1, 2, M (3-14) 

m-1 P" 1 m p 

where M Is the total number of grid points. The solution of the p equation 
with respect to the central unknown 4 p yields: 



This equality will not be satisfied until the problem has relaxed or converged 
to the final result. When (3-15) Is not satisfied, we assume the rlghthand 
side to be the better value for so by denoting the various approximations 
by the Index o we step through the Index p, replacing with ap 0 * 1 to arrive 
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at the computational rule* 






M 

I 


n*p+l 




We have jumped ahead and added the over-relaxation factor W. In the ordinary 

single step method* W ■ 1. and for W > 2 the method will diverge. The matrix 

c , obtained from the discretization of a boundary value problem Is of the 
pm 

banded symmetric-definite type* and for such* convergence Is Insured for 
W < 2; Panel uses W ■ 1.9 with no divergence problems. 

In the program, the subroutine FIELD controls the Poisson calculation. 
The calculation of the Interlur coefficients 1$ delegated to the subroutines 
CNS, uEW* and CUD. The boundary conditions and B.C. Influenced coefficients 
are effected In FIELD and the subroutine RELAX performs the relaxation 
operation. 


THE VLASOV PROBLEM 

In kinetic theory, the density and current at a point x' are given by the 
Oth and 1st velocity moment of the single particle distribution function: 


N x (x* ) ■ Jfj(x', v* ) d 3 v‘ 

(3-16) 

l f (*> ■ % *') V'*n d 3 v\ 

(3-17) 


where n Is the unit vector In the direction of The distribution function 
Is the density of particles In six-dimensional phase space (three position and 
three velocity coordinates). Further progress now requires flndlm. ' at x' . 



Application of Llouvllle's theorem to a colllslonless plasma leads the col 11- 
slonless Boltzmann or Vlasov 1, equation (Montgomery and Tldman, 1964). 

■ £ ♦ ».vf, ♦ (f ♦ . 0. (3-18) 

In words, f, Is constant along a particle's path In six-dimensional phase 
space, which can be characterized by the constants of the motion. In a 
general electrostatic field such a constant Is the total energy, defined by 

H,(x, v) ■ £ m s v 2 ♦ q s V(x) 

where q $ V(x) Is the potential energy of the particle at x. The six-dimension 
phase space path projected onto the usual space coordinates Is just the usual 
trajectory of a particle prescribed by Newtonian mechanics. 

Consider the trajectory connecting (x', v * ) with (x, v) for a given elec- 
trostatic field where at x, the distribution of particles of specie s Is known 
to be f s (x, v). If fj can be written as a function of only H(x, v), and since 
H(x' , v') • H(x, v) ; we have therefore 

f s (H(x, v)) ■ fj(H(x\ v')); 


t The Vlasov equation represents the zeroth order terms In a cluster 
expansion of the Llouvllle equation, with smallness parameter g ■ (np^j 3 )" 1 , 
the Inverse of the number of particles In a Debye sohere. For qEO under 
substorm conditions n e • 1/cc and kT. • lOkev, g ■ 10~*, so the colllslonless 
approximation Is a very good onr.. However, In the F region with n • 10 6 /cc 
and kT e ■ .2ev, g - 10 suggesting that the transport problem there needs a 
more detailed treatment. 


or simply 


(3-19) 


f;(x\ v') ■ f,(x, v). 

We have now effectively solved the Vlasov equation and may now, In principle, 

3-ifc 3-<7 

proceed to evaluate the Integrals In equations 4 and 9 , 

Note that some different value of v' will map to a different point 
(x 2 . v 2 ) where we know the distribution function to be different (or zero). 
Thus, In evaluating the Integrals In equations (3-16) and (3-17), equation 
(3-19) must be used to develope a composite expression for f. For example, 
consider the problem of a non-emitting body immersed In a Maxwellian plasma. 
At Infinity, the distribution fun-cion In three-dimensions, for specie s Is, 

fj(*. v) ■ Ns( 2 Jfc ) 3/2 «xp [- (3-20) 

At some point x' near the body, the distribution function will be, 

f$(x\ v') ■ exp [-(y mjVj 2 ♦ q,V(x'))AT s ] * G,(x‘, v'), (3-21) 

where G f Is a function with a value of either zero or one depending on tdiether 
(x‘, v 1 ) maps to a non-source or source at Infinity. In other formulations, 
the G function Is effectively replaced by reconstructing the limits of 
integration In equations (3-16) and (3-17). 

In practice, the Integrals In (3-16) and (3-17) are approximated by 
summations over a discrete set of velocities where each value of v* represents 
a trajectory that must be followed to evaluate G(x', v'). We now have the 
choice of either starting trajectories at "Infinity" and following them In; or 
because of the assumed time-independence, we could start at x' and follow 
trajectories backwards In time to "Infinity". The first technique has been 
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dubbed the "Outside-In Method" by Parker (1964) and has the advantage of 
having all trajectories successfully connecting to a source and of supplying 
useful 1 trajectory Information to all points along the trajectory. Its chief 
disadvantage lies In the difficulty of getting adequate trajectory probing of 
some regions of the problem. The Inside-out Method adopts the other approach 
of following trajectories backwards In time. This allows one to evaluate 
G(x' , v 1 ) at all points with equal accuracy, but can lead to large numbers of 
trajectories to be retraced with each Iteration. This last difficulty has 
been recently overcome by recording the fate of each trajectory so that In 
subsequent Iterations, that Information can be used to trace only those 
trajectories that lie on the velocity space boundary between null and escaping 
trajectories. This "boundary tracking" Innovation can greatly Increase storage 
requirements, but the reduction In time requirements make It essential. 

In the following paragraphs I shall described how PANEL performs the 
Integrals (3-16) and (3*17). Parts of this description has been taken 

directly from Parker (1977). 

It Is convenient to transform (3-16) and (3-17) to energy and angle 
variables In velocity space. Since we will be primarily Interested In 

Maxwellian energy distributions, we may adopt the following units In terms of 
which dimensionless variables may be defined: 

kT ■ unit of energy, whera ~ Is the temperature of the Maxwellian 
distribution; 

/2KT7m ■ unit of velocity, namely, the most probably thermal velocity; 

Nq ■ unit of particle density, the unperturbed density; 

J Q • fTTJ7* m ■ unit of current, the undisturbed thermal current. 

The energy and angle variables are: 

H ■ energy In multiples of kT; 
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a ■ polar angle with resptet to z-axls; 

$ ■ az<«utha1 angle with respect to tha plana containing tha z-axls and 
tha point x. 

Thasa angles which define tha orientation of tha velocity vector v, are 
Illustrated In figure 3/2. Note that tha potential energy a Is also In units 
of kT. so with tha new unit of velocity we can write: H ■ v* ♦ t. 

The density and current Integrals (3*16) and (3-17) may be written 

N # ■ /// f'v*dv sin a da dp G 

J $ ■ /// f'v*dv cos a slna da dp G 

where J s Is assumed to be the current to a surface perpendicular to the 
z-axls. Introducing now. the Maxwell distribution (without drift), we have 

n - it - -4 tt r dH f* Sin a da J 2t G dP (3-22) 

1 • ^ ■ i- r t’ H (H - •) dH J* /2 cos a sin a da J 2w G dp (3-23) 

w aax(O.e) 0 0 

The lower Halt on the energy Integral Is chosen to be zero for an attractive 
potential (♦ < 0), and to be a for a repulsive potential (a > 0). This 

ensures that we never consider particles with negative kinetic energy. 

The Integrals In (3-22) and (3-23) are evaluated by the Method of 

Gaussian quadrature (Jennings. 1964). It Is not feasable to derive It here, 
but the wthod can be Illustrated by stating that In the foraula 

/J f(x) dx ■ l A.ffx.) + R 
1-0 
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It Is possible to choose and x 1 such that R ■ 0 for f(x) any polynomial of 
degree < 2n + 1. When attempting to Integrate a function of undetermined 
degree, It Is desirable to make n as large as possible. This poses two 
problems: 1) the formulas for x^ and Aj get Increasingly complicated as n 
Increases and 2) the step function nature of G(x, v) Implies a polynomial of 
Infinite degree. Both problems are partially overcome by dividing the 
Integration Interval uniformly as for ordinary trapezoidal Integration, and 
then applying a Gaussian quadrature of order 2 to eac! subinterval. Thus, In 
preparation for this, we transform the ranges of Integration Into Intervals 
between -1 and +1 by the transformations: 


n(c) *747* (0. ♦). 


-1 < e < *1 


(3-24) 


a * COS' 


a * cos -1 (-a) 


for current 
for density 


-1 < a < +1 


(3-25) 


6 * ir (1 + b). 


-1 < b < +1 


(3-26) 


The transformed current and density Integrals then become. 


„ . J_ ; + i j+i e - H(c) / H[c) - * g da db dc . 

n 1 1 1 (i - c) 2 


(3-27) 


j • 7 /!} /!} /!} e- H(c) (H(c) - ♦) 


G da db dc 
(1 - c) 2 


(3-28) 


We now have the Integrals In a form suitable for Gaussian quadratures. 
We now divide the c-range Into M c sub-intervals, and apply a Gaussian 
quadrature of order 2 to each. Similarly, we divide the a-range and b-range 
Into M a and sub-intervals, with Gaussian quadratures of order 2 applied to 
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each sub-interval. Now, both equations (3-27) and (3-28) can be put In the 
form 


1 ■ /!} £} /!}«(H)"6(H. a, 0)*da db dc 


(3-29) 


which may be approximated by the sum: 


Mj^Mc 


u |y| |yt 

f [ a ? [W(H>G(H*. a*. 0") ♦ W(H + ).G(H*, a + , 0 + )] 

Kc-1 Ka*l Kb*l 


(3-30) 


where W Is the energy weight function defined by. 


W(H) 


e' H(c) r H(c> - *1 


W(H) 


2(1 - c) 
e~ H(c) 


for current 


(3-40) 


/rT (1 - c)" 


for density 


with 


H" » H(c“) , H + = H(c + ) 
a" « a(a“) , o + * a(a + ) 
6* * 0(b") , 0 + * 0(b + ). 


Gaussian quadrature of order 2 applied to the Interval -1, +1 yields the 
— 1 /? 

abscissas i (3) with a weight coefficient of unity. Applying this simple 
formula to each sub-interval gives the formula, 


i* ’ *7 [* 7? + 2K i ‘ i * M, 1 for 1 ‘ *• b> c - 


(3-41) 


so, with these formulae for a, b, and c, we may use equations (3-24) through 
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\ 


(3-26) to choose sets of trajectories to be followed backwards In time to 
either source or non-source regions and thus approximate the Integrals for 
density and current. 

In the first Iteration only a small set of the total 8M c M a M b trajectories 
are used, and In following Iterations the number Is Increased until the 
maximum number of trajectories Is being followed. This Is an economizing move 
that doesn't affect the accuracy of the calculation since In the beginning 
Iterations, densities are only approximate. As each trajectory Is followed to 
Its end-point, Its fate (escape ■ true, absorbed ® false) is recorded In the 
four-dimensional logical matrix (N x 2M C x 2M a x 2M b ) called TRYE for elec- 
trons and TRY I for ions, where N identifies the point. When each trajectory 
has been used at least once, the TRY matrices are used in subsequent Itera- 
tions to trace only those that lie on a velocity space boundary. This is 
accomplished by simply comparing the last recorded fate for the trajectory in 
question to that for each of its six energy and angle neighbors In the TRY 
matrix. If all seven fates are Identical, the trajectory is not followed and 
its fate is merely read from TRY. If they are not all Identical, the trajec- 
tory is followed and any change of fate Is recorded In TRY. If all the fates 
for a point and particle become the same, a shut-out would occur and no 
trajectories would be traced. To prevent this, all of the highest energy 
trajectories are exempted from the boundary searching process, and traced each 
Iteration. 


i 


i 

i 
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Trajectories 

In an electrostatic field with no magnetic fields, a particle will move 
according to the equation. 


X + Jt +4£ft 2 

o i m 


where t Is time and £ Is the electric field. This equation Is put Into the 
units common to PANEL with the transformation t ♦ t ' • (2kT/m) . Thus we have 


x o + V + \ 


*0 * V* + 1 


z o * V + \ 

K-7& t<2 


(3-42) 


where as before, v Is In units of the most probable thermal velocity, and * Is 
the dimensionless potential. 

PANEL traces particle trajectories on the same grid that Is used for the 
Poisson calculation. At each Interior grid point, the six neighboring Inter- 
mediate points each define a face of a cell enclosing the grid point. The 
velocities In equations (3-42) are always given as a result of a previous 
step, or as initial conditions as a trajectory starts. The partial deriva- 
tives of the potential In (2-42) are approximated by divided differences 
calculated In the following manner. At the point P(x^, yj, z^), form the west 
and east potential differences. 


4 *W * *(1. J, k) • j, k) 


4 *E ' j, k) ‘ *(1. j. k) 
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and similarly for y and z. Next, the absolute value | (a* w - a« e )| Is compared 
to the particle's total energy, H, multiplied by an Input resolution factor 
RES. If the variation In potential differences Is less than H • RES, PANEL 


uses 


AX " X1 +1 - XI. 1* 


for the entire cell, and If the variation Is too great, the cell Is halved and 
for the west and east halves we use. 



A* 


X - X 
1 1-1 


and 



Thus, a cell can conceivably be divided Into eight sub-cells. This subdi- 
vision is always performed at the start of a trajectory when the particle will 
probably never enter most of the sub-cells. 

Once the cell or sub-cell has been defined and the "electric field" 
calculated equations (3-42) are solved Independently for the times required to 
cross the cell, and the shortest positive time Is chosen. Using this time, 
the particle Is stepped to another face of the cell and the process begins 
again In the next cell until an outer boundary Is reached. At a boundary, a 
particle can escape, be absorbed on a surface, or be reflected (for a 
reflection boundary condition). 
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STABILITY AND CONVERGENCE 

Parker and Sullivan (1970) have anallzed the stability and convergence of 
the Inside-out method, applied to a uniformly charged sphere In a uniform 
plasma. Although a three-dimension method like PANEL could be expected to 
differ from a simple one-dimension method In Its stability properties, tests 
have shown that the results of their study are applicable to PANEL. That 
analysis will be briefly outlined. 

IT we Imagine that the Poisson solving process can be represented by the 
operator L(*) and that the Vlasov process can be represented by £(*), then the 
state of a system would be prescribed by 

L(«) ■ F(#), (3-43) 

and the Iterative procedure previously described would follow the Picard 
iteration rule, 

L(* n+1 ) * F (* n )» (3 " 44) 

where n Is the iteration Index. This Iteration scheme, however, was found to 
diverge when the distance between their sphere and the model boundary exceeded 
the Debye length. An effective cure for this, Is to replace rule (3-44) with 

L(Vl> ' F (“*n + (1 * «) O- 0 < . < 1. (3-45) 

This technique Is called mixing, and the superscript M indicates previously 
mixed potentials, l.e.. 
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ip * 1 m 


♦U-l • ®*n-l + » * *)«n-2 


Their analysis of rule (3-45) predicts monotone convergence for 
a < 2/(2 ♦ y), oscillatory convergence for 2/(2 + y) < a < 2/(1 ♦ y), and 
divergence for o > 2/(1 + y), with an optimum value, ■ 2/(1 ♦ y). The 
parameter y Is given by 

y ■ 2d z /»% 2 , 

S 

i 

! 

J 

where d Is the boundary-object separation distance, and Is the Debye length 
for the plasma. These predictions have been proven to be accurate for the 
one-dimension sphere model. 

One could Insure convergence by choosing a very small, but then a large 
number of Iterations would be required. Therefore It Is deslreable to 
optimize o. With PANEL, I have found this prediction for to be a good 
first approximation; but to really obtain optimal convergence, calculations 
must be stopped every three or four Iterations to adjust a. 
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4: RESULTS 


The results presented here are of two distinct types: "calibration" 
models designed to test PANEL against problems for trfilch analytic answers are 
available* and two production models to demonstrate PANEL'S capabilities. 
Thus far, the emphasis has been on the former. Just as with an Instrument, 
results are worthless without calibration. This emphasis has been rewarded, 
as many subtle errors (both with PANEL and ny use of PANEL) have been detected 
and corrected. For this reason, the runs presented here represent only a 
fraction of those that have been made. 

The models called Gap 06, Gap 07, and Gap 08 are calibration models of 
the problems described In equations (2-12) through (2-17) In section 2 of this 
report. Pan 21 Is a model of a planar electron diode, and can be compared to 
the Chlld-Langmulr law, equation (2-20). Finally, Pan 29 and Pan 36 are 
production models of charged panel In a plasma similar to that encountered In 
the Chamber A experiments at the Johnson Space Center (McCoy and Konradl, 
1978). All of these are two-dimensional models. Three-dimensional tests have 
also been made, but limited computing time has preveri.^ the running of 
physically meaningful three-dimensional models. 

In the gap problem, electrons are accelerated from a cold cathode (T • 0) 
to the potential Vi (see figure 2/1) at x ■ 0, to produce a beam current J. 
PANEL models this experiment by assuming that there Is an undisturbed 
Maxwellian plasma at < -d, so that the current J Is the randor. thermal 
current (J Q * N o e /kf/Zvm) crossing the grid at x - -d. (Although the plasma 
has density the electrons crossing the grid have N ■ Nq/2). As described 
In section 2, this current Is normalized by the relevent C-L current, thus we 
have the current ratio. 
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0 ■ J/J CL « 

The results of Gap 06, Gap 07, and Gap 08 are plotted In figures 4/1 and 
4/2. In these plots, the transmitted electrons travel from right to left 
across a gap of one meter. This Is modeled by 24 grid points; 12 z and 2 x 
coordinates. At z - 0 electrons are absorbed; at z ■ 11 (/iot shown), they are 
generated; and they are reflected at both x boundaries. (Since this Is a one- 
dimension problem, PANEL could have been fitted with a one-dimension option, 
but unlike the two-dimension option, a one-dimension option would have only 
limited applications.) In all three plots, the potentials predicted by the 
classical theory are labeled as curve A, and the results of PANEL are labeled 
P. The features of these models are: 

Gap 06: 0 » 10, J - 2.373 x 10‘ 2 A/m 2 , T 0 - 10 eV. Vi - 100 V. 

N 0 » 2.8 x 10 5 cm' 3 , M 0 - 4, M a - 32 ; 

Gap 07: 0 - 10, J • 2.373 x 10" 2 A/m 2 , T e • 1 eV, Vi - 100 V, 

N q - 8.9 x 10 5 cm' 3 , M g « 4, - 32; 

Gap 08: 0 ■ 4, J ■ 9.49 x 10' 3 A/m 2 , T e » 1 eV. - 100 V, 

N q - 3.5 x 10 5 cm' 3 , M e ■ 4. M a - 32. 

Gap 07 and 08 are both well converged, but Gap 07 has an uncertainty Indicated 

by the error bar on the plot. I consider these models to be a positive test 

of PANEL, inspite of the large deviations from the classical predictions. The 
classical theory considers a source of electrons with no thermal spread. By 
comparing Gap 06 with Gap 07 we can see that as the source plasma cools from a 
temperature of 10 eV to 1 eV, the results get closer to the classical predlc- 
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tlon. In Gap 08 where 8 ■ 4, the predicted minimum potential Is 75.0 volts 
while PANEL gives 75.1 + 0.7 (the error Indicates the degree of convergence). 
This again Indicates that the disagreement with the classical theory In Gap 06 
and Gap 07 are due to non-zero temperatures since one would expect this effect 
to be most pronounced with low minimum potentials, and least pronounced with 
higher minimum potentials. 

From Gap 08 we can also learn something about the number of trajectories 
that must be traced to give accurate densities. In figure 4/2, the lower 
curves labeled D and C are densities for PANEL and classical theories respec- 
tively. Briefly, the classical densities are derived by eliminating v from 

J • Nev (4-1) 
and j mv 2 ■ eV (4-2) 
to get N • J( m /2e 3 V) l/2 . (4-3) 

C 

For Gap 08 the total zenith angle range of 2w Is covered by 64 trajectories to 
give a trajectory separation of .098 rad. or 5.63°. This separation was 
further reduced by one half by noting that due to symmetry, positive and 
negative angles of equal magnitude lead to equivalent trajectories; thus all 
trajectories were shifted by half of the separation angle. The result Is that 
although the voltages were obtained with good accuracy, the densities still 
lack resolution. 

Pan 21 represents a simple but Important test of PANEL. This Is a 
comparison of PANEL with the Chlld-Langmulr law shown In fig. 4/3. Due to the 
close agreement, a curve has been drawn only through the PANEL points. At 
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selected points, PANEL and C-L potentials are given for comparison. The C-L 
.o.entlals are given In parentheses and C-L densities are plotted with 
crosses. Here 32 points (2 x 16) were used to model a diode with a 16.51 
meter plate separation, and a 100 volt potential difference. The model 
parameters are: 

! 

Pan 21; T # ■ 1 eV, N ■ 3.2 x 10 2 cm* 3 , X d ■ 0.4m, \ ■ 4, M fi ■ 32, 

J- 8.58 x 10“ § A/m ? 

The greatest disagreement between PAN 21 and the C-L theory occurs at z ■ 

14, whete the PANEL prediction Is 22t high, with Improved agreement at lower z 
values. At z ■ 8, the disagreement Is only IX. The larger deviations should 
be expected In the low voltage region near the cathode due to the non-zero 
Injection velocity of the electrons. For this reason. It would be deslreable 
to compare PANEL predictions with the modified C-L law, equation (2-10), but 
unfortunately, this comparison has not yet been made. 

Pan 29 and Pan 36 are two-dimensional models of a cross-section of an 
Infinitely long, one meter wide panel held at a potential of 100 volts In a 
hydrogen plasma with equal Ion and electron temperatures of 10 eV. The chosen 
plasma temperature of 10 eV Is higher than the usual temperatures encountered 
In LEO or In the JSC Chamber A experiments which are frequently less than 1 
eV. Models with a panel potential of 100 V and a temperature of 1 eV have 

been considered, but under these conditions PANEL Is significantly less 
stable. To achieve stability thus requires a smaller mixing parameter, more 
Iterations and more computing time; so for these first models, a higher but 
not unreasonable temperature was chosen. 

Due to the symmetry of the problem, It was possible to model the entire 
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cross-section by calculating potentials and densities In one quadrant only by 
using the reflection boundary condition on the DOWN and WEST boundaries. For 
PAN 29 the UP B.C. is V ■ 0, and for the EAST boundary the B.C. Is the zero 
normal gradient boundary condition [B.C. (2), iV/an ■ 0]. Model parameters 
for PAN 29 and PAN 36 are: 

T 1 • 10 eV, T # - 10 eV, • 1.9 x 10* cm’ 3 , \ Q ■ .17m, 

1.61 x 10’ 3 A/m 2 , ME - 4, MA - 32. 

PAN 29 potential contours are displayed In Figures 4/4, and the PAN 36 results 
are displayed In Figures 4/5, 4/6 and 4/7. The potential contours shown In 
Figures 4/4 and 4/5 were produced by the well-known technique of eyeball 
Interpolation. For PAN 36, potentials (labeled P) and charge densities 
(labeled C) along the DOWN and WEST boundaries are presented In Figures 4/6a 
and b respectively, and the electron and proton densities along these same 
boundaries are shown In Figure 4/7. In addition to the change In B.C.'s the 
other differences between PAN 29 and PAN 36 are location of the boundaries and 
the number of grid points. PAN 29 has an UP boundary distance of 1.85 meters 
with a 8 x 8 grid, while PAN 36 has an UP boundary distance of 2.4 meters with 
a 9 x 10 grid. 

Several Interesting features of the general problem can be observed by 
comparing the contour plots 4/4 and 4/5. First, by assuming that PAN 36 Is 
the better model of the two, we can see that B.C. (2) on the EAST boundary of 
PAN 29 (Figure 4/4) produced potentials along the DOWN boundary that are very 
close to those of PAN 36; and therefore supposedly better than tftat could have 
been obtained with a V ■ 0 condition. The drawback Is that B.C. (2) severely 
destabilizes the problem. One model, PAN 35 (nearly Identical to PAN 36, but 
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with B.C. (2) on both the UP and EAST boundaries) diverged severely with a 
mixing factor of 0.08, while PAN 38 was stable with mixing factors up to 
0.5. This suggests the possibility of speeding convergence on smaller grids 
by using varying mixing factors with small values near B.C. (2) boundaries and 
Increasing to larger values near fixed potentials. 

It can also be seen that the Increased grid point density near the panel 
In PAN 36 has caused the 70 V, 50 V, and 30 V contours to move closer to the 
panel with smoother contours near the edge of the panel. The electron 
currents collected from above the panel are indicated by the arrows below the 
panel In figures 4/4 and 4/5, and have been normalized by the random thermal 
current* J Q . Near the center of the panel* both models give the same current* 
but with the Increased point density In PAN 36 we begin to see a slight 
reduction In current collection near the edge. A further Increase In point 
density would probably show more current focusing. However* the strong 
central focusing (greater than an order of magnitude difference between 
central and edge currents) observed by McCay (1980) In the solar panel tests 
at JSC Is not Indicated In these models. This focusing could be dependent on 
the correct choice of panel voltage and plasma parameters* but Is most likely 
due to a band of dielectric along the edges of that test panel. This 

possibility will be tested In future models. 

For both PAN 29 and PAN 36* the C-l screening distance Is D CI _ - 1.2 m, 
and the corrected screening distance Is D s ■ 1.73 m. These points are 
Indicated In Figures 4/4* 4/5* 4/6, and 4/7 . and the uncorrected C-L contour 
Is marked with crosses In Figure 4/6. Figures 4/5 and 4/7 show that poten- 
tials have been reduced to less than kT/e (■ 10 V) within either estimate. 
There Is some "compression 11 of the contours caused by the closeness of the 
V ■ 0 boundaries, as Is evidenced by the most distant points In Figure 4/7 
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where the electron density unrealistically drops below the proton density due 
to the artificially high electric field between the outermost two points. 
However* comparison of the EAST boundaries of Figures 4/4 and 4/S suggests 
that this concession Is not too severe. 

Although all of the PAN 29 and PAN 36 boundaries are too close to allow 
an undisturbed plasma region to develop* Figure 4/7b shows a definite pre- 
sheath region beyond the point with electron and oroton densities nearly 
equal but reduced from the ambient values. * 

The Models PAN 29 and PAN 36 are clearly not a complete study of the 
solar panel -plasma interaction problem, but the results that have been 
presented should demonstrate that PANEL Is capable of accurate space charge 
calculations. Three-dimensional test calculations have been successfully run, 
but time limitations have prevented full scale three-dimensional modeling. 
The two-dimensional model PAN 29 required twelve minutes of processing on an 
ITEL AS/6. PAN 36 required about thirty minutes, and three-dlmenslonai models 
are expected to require many hours of processing time. Although this Is not 
an extreme requirement, two-dimensional modeling will continue to be Important 
for deciding matters such as the placement of boundaries, or the effect of 
chamber walls. 
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Appendix A: Symbols and Constants 

The sy.uools and constants used throughout this report are reviewed In 
this appendix. With only a few noted exceptions, the M(S system of units has 

been adhered to. 

I 

T; temperature : i degrees Kelvin. 

k; Boltzman's constant » 8.62 x 10’ 5 eV/°K. 

e; electon charge ■ 1.602 x 10’ 19 Coulomb. 

m e ; electron mass » 9.11 x 10’ 31 kilograms. 

-27 

m p ; proton mass * 1.673 x 10 grams. 

e Q ; free space permittivity » 8.85 x 10’ 12 farad/meter. 

V; potential in Volts. 

dimensionless potential normalized by kT/e. 
v; velocity In meters/second, or normalized by /2lcT/m. 

E; dimensionless electric field. 

H; total energy In Joules, or dimensionless total energy normalized by 
kT/e. 

N; density in m’ 3 . 

N Q ; ambient density in m -3 . 
n; dimensionless density normalized by N Q . 

J; current In Amperes/meters * Am . 

J 0 ; random thermal current in Am . 
j; dimensionless current normalized by J Q . 
p; charge density In Coulomb m’ 3 . 
a; sheath thickness. 
x Q ; Debye length. 
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0RBIT3: Individual 
trajectory incre- 
ments 


INTERP: Locate 
particles after 
trajectory steps, 
calculate E fields 




















Appendix C: The 2-D Option 


When a problem exhibits sufficient symmetry. It Is sometimes possible to 
reduce the number of Integrals In equations (3.16) and (3.17) that must be 
performed numerically; so. In a sense, PANEL'S 2-0 option still produces 3-0 
results. 

Using the Maxwellian distribution given In (3.21), and writing (3.16) and 
(3.17) in terms of Cartesian velocity coordinates we have, 

N < x '> * N 0 fe )3/2 5 v i (c-u 

« -p t- m ♦ *; 2 ♦ *;*) - 

3 (x-) - q n o C^ kt J 3/2 £; /3v; /3v* g(x\ v) (c-2) 

* (1 • ») e*P 1- (v; 2 + v^ 2 + vj 2 ) - 

If G(x', v') is not a function of v y , and if n in (C-2) lies entirely in the 
x-z plane, the Vy integrations can be performed immediately leaving 

n(x') * ViSet) &z g( x‘' v> (C-3) 

■ «P [ - * K* . v; 2 ) - 

and 

J z (*■) ■ N o^) /!*; />; G < x '- v') v z (c-4) 

«p t- * w; 2 ♦ *; 2 ) - 
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where we have set n ■ z. Transforming now to cylindrical polar coordinates, 
and as In Section 3, denoting the exponential argument by H; we have 

N(x') * N 0 & /^v'dv'J** de ' G(x ' , v') exp (-H) (C-5) 

and 

J z * N o *‘Vj3 s e' de' G(x' , v* ) exp (-H) (C-6) 

Again as In Section 3, we transform to dimensionless v, H, and *, and norma- 
lize by Njq and J so to get the two-dimensional versions of (3-22) and (3-23); 

n * 4r /‘e- H dH / + > G(x' , J' ) (C-7) 

* max(0,*) 

and j 3 — - J" e" H /H - * dH J ^cose de G(x' , v' ) (C-8) 

z max(0,*) - 11 / 2 

In preparation for Gaussian quadrature, we make the following transformatins : 


E = | + max(0,*). 


-1 < c < +1 (C-9) 


9 * ira 

0 3 s1n” l (a) 


,for density 

-1 < a < +1 

,for current 


(C-10) 


Including these we have, 


and 


n * J*} r} e~H(c) <*= da G (x‘, »') 

** (1 - c) 2 


j * — /+} /♦} e- H < c > /irrr^^Gtx'. ?•> 
2 G _1 (1 - c) 2 


(C-ll) 


And, as in section (3), we Introduce the sub-interval Gaussian quadrature 
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approximations: 


M M 

" or j 2 * ffV £ C ^ [M 2 (H-).G 2 (H-..-) ♦ W-(H + )-G,(H + ,0 + ) (C-12) 

z c "a K «1 K »1 22 22 

c a 

where the two-dimensional energy weight function Is 

H , e* H ( c ) x | 1 ,for density 

2 (1 - c)2 /h - $ ,for current (C-13) 

/ * 

and H* » H(c’), e i » 9(a i ). 

The Gaussian abscissa formula (3-41) and the transformations (C-9) and (C-10) 
are used to initiate trajectories. 
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FIG. 2/1 

Schematic representation of a planar electron 
diode, showing the effect of the electron space charge 
on the interelectrode potential for two different 
arrangements of the electrode potentials. 




Potantiol </4C> 


Input courrent density 6 

Fig.2/3a. Minimum potential aa a function of input 
current 3* The dotted portion ia called the "C-overlap H 
solution and i8 excluded by atability considerations. 



Fig.2/3b. Potential profiles 
in the interelectrode space for 
the classical, positive potential 
minimum solutions with input 
currents between 0 and 8 . 
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Fig.2/3c. Potential profiles 
in the interelectrode space for 
the classical virtual-cathode 
solutions ( ♦a ■ 0 ) with input 
currents between 4 and « . 






FIG. 3/1. A sample grid illustrating a projection 
onto the X-Y plane of the volume associated with 
the point P, and the directions indicated by N,S, 

E § W t U t and D . 



FIG. 3/2. An illustration of the polar- spherical 
angles a and 0 used in PANEL. 
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FIG. 4/1: In both graph* , curve P ia PANEL'S prediction for 
the inter-electrode potentials/ curve A is the classical 
prediction/ and curve D is PANEL'S prediction for the 
electron density. The Z distance unit is 0.1 meter. 

6 is the normalized current density. 
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FIG.4/2: The upper curve is PANEL'S prediction for the 
inter-electrode potentials, and the classical potentials 
are the unconnected crosses. Of the lovrer curves, 0 is 
the PANEL result for the densities, and C gives the 
classical densities. The Z unit of distance is 0.1 meter, 
and B is the normalized current density. 
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FIG. 4/3: PANEL'S predictions for the Child- Langmuir 

electron diode are plotted with the connected dots, 

and selected values of potential and density are 

presented. The classical C-L density predictions are 

plotted by the unconnected crosses, and for comparison, 

selected values of the C-L potential and density are 

given in parenthesis. For this model, N o * 320/cc, 

- 6 2 

J 0 ■ 8.6x10 A/m , and the total diode separation is 


16.51 meter. 
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FIG. 4/4, Model Pan29 : Equipotential contours about aim 
wide, infinitely long panel, uniformly charged to 100 V ; 
modeled with 64 grid points (not shown). The plasma 
temperature and density are kT 3 10 eV , N Q * 1.9 x 10^/cc. 

The thermal currents are J Qe ■ 1.61 x lO'^A/m^ for electrons, 
and J Qp ■ 3. 76 x 10" 5 A/ra 2 for protons. The electron currents 
collected by the panel are indicated by I’he arrows below 
the panel and are normalized to J ; the average proton 

_q oe y 

current to the panel is 1.5x10 A/m . The corrected and 
uncorrected Child-Langmuir screening distances are indicated 
above and are respectively: D s 3 1.73m, D CL « 1.2 m . The 
X and Z unit of distance is 1.0 meter. 




FIG. 4/5 , Model Pan36: The same panel and plasma of model 
Pan29 are modeled here with 90 grid points (not shown) , 
and with a V = 0 boundary condition on both the UP and 
EAST boundaries. The distance unit is 1.0 meter. 
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FIG. 4/6a 



FIG. 4/6: Potential (curve P) and charge density (curve C) profiles 
alonge the DOWN (4/6a) and WEST (4/6b) boundaries for the model 
Pan36. The density scale on the right is normalized to the ambient 
density, N 0 ■ 1.9 x loVcc. In 4/6b, the Child-Langmuir potential 
profile is Indicated by the unconnected crosses, and the initial 
velocity corrected C-L screening distance, D fl , is indicated by 
the arrow. The distance unit is 1.0 meter. 
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FIG. 4/7: Potential (curve P) , and electron and proton density 
profiles along the DOWN (4/7a) and WEST (4/7b) boundaries for the 
model Pan36. The corrected and uncorrected Child-Langmuir screening 
length estimates D s and D^l are indicated in 4/7b. The distance 
unit is 1.0 meter. 
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